!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2015 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!     
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of 
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine objective_shapeener_dx(co,kon,ipkon,lakon,ne,
     &  stx,elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,
     &  ielmat,ielorien,norien,orab,ntmat_,t0,t1,ithermal,prestr,
     &  iprestr,iperturb,fn,iout,qa,vold,nmethod,
     &  veold,dtime,time,ttime,plicon,nplicon,plkcon,nplkcon,
     &  xstateini,xstiff,xstate,npmat_,matname,mi,ielas,icmd,
     &  ncmat_,nstate_,stiini,vini,ener,enerini,istep,iinc,
     &  springarea,reltime,calcul_qa,iener,ikin,nal,ne0,thicke,
     &  emeini,pslavsurf,pmastsurf,mortar,clearini,nea,neb,ielprop,
     &  prop,distmin,ndesi,nodedesi,ndirdesi,nobject,g0,dgdx,
     &  numobject,sti)
!
!     calculates stresses and the material tangent at the integration
!     points and the internal forces at the nodes
!
      implicit none
!
      integer cauchy
!
      character*8 lakon(*),lakonl
      character*80 amat,matname(*)
!
      integer kon(*),konl(26),nea,neb,mi(*),mint2d,nopes,
     &  nelcon(2,*),nrhcon(*),nalcon(2,*),ielmat(mi(3),*),
     &  ielorien(mi(3),*),ntmat_,ipkon(*),ne0,iflag,null,
     &  istep,iinc,mt,ne,mattyp,ithermal(2),iprestr,i,ii,j,k,m1,m2,jj,
     &  i1,m3,kk,iener,indexe,nope,norien,iperturb(*),iout,
     &  nal,icmd,ihyper,nmethod,kode,imat,mint3d,iorien,ielas,
     &  istiff,ncmat_,nstate_,ikin,ilayer,nlayer,ki,kl,ielprop(*),
     &  nplicon(0:ntmat_,*),nplkcon(0:ntmat_,*),npmat_,calcul_qa,
     &  nopered,mortar,jfaces,iloc,igauss,desvar,ndesi,nodedesi(*),
     &  ndirdesi(*),nobject,numobject 
!
      real*8 co(3,*),shp(4,26),stiini(6,mi(1),*),
     &  stx(6,mi(1),*),xl(3,26),vl(0:mi(2),26),stre(6),prop(*),
     &  elcon(0:ncmat_,ntmat_,*),rhcon(0:1,ntmat_,*),xs2(3,7),
     &  alcon(0:6,ntmat_,*),vini(0:mi(2),*),thickness,
     &  alzero(*),orab(7,*),elas(21),rho,fn(0:mi(2),*),
     &  fnl(3,10),beta(6),q(0:mi(2),26),xl2(3,8),
     &  vkl(0:3,3),t0(*),t1(*),prestr(6,mi(1),*),
     &  ckl(3,3),vold(0:mi(2),*),eloc(9),veold(0:mi(2),*),
     &  springarea(2,*),elconloc(21),eth(6),xkl(3,3),voldl(0:mi(2),26),
     &  xikl(3,3),ener(mi(1),*),emec(6),enerini(mi(1),*),
     &  emec0(6),veoldl(0:mi(2),26),xsj2(3),shp2(7,8),
     &  e,un,al,um,am1,xi,et,ze,tt,exx,eyy,ezz,exy,exz,eyz,
     &  xsj,qa(3),vj,t0l,t1l,dtime,weight,pgauss(3),vij,time,ttime,
     &  plicon(0:2*npmat_,ntmat_,*),plkcon(0:2*npmat_,ntmat_,*),
     &  xstiff(27,mi(1),*),xstate(nstate_,mi(1),*),plconloc(802),
     &  vokl(3,3),xstateini(nstate_,mi(1),*),vikl(3,3),
     &  gs(8,4),a,reltime,tlayer(4),dlayer(4),xlayer(mi(3),4),
     &  thicke(mi(3),*),emeini(6,mi(1),*),clearini(3,9,*),
     &  pslavsurf(3,*),pmastsurf(6,*),distmin,enerint,g0(numobject),
     &  dg,dgdx(ndesi,nobject),sti(6,mi(1),*)
!
      include "gauss.f"
!
      iflag=3
      null=0
!
      mt=mi(2)+1
      nal=0
!
!     -------------------------------------------------------------
!     Initialisation of the loop for the variation of 
!     the internal forces
!     -------------------------------------------------------------
!
!     Loop over all designvariables
      do desvar=0,ndesi
!
!        initialisation of dg
! 
         dg=0.d0
!
!     Loop over all elements in thread
      do i=nea,neb
!
!        initialisation of enerint and 
!
         enerint=0.d0
         lakonl=lakon(i)
!
         if(ipkon(i).lt.0) cycle
!
!        no 3D-fluid elements
!
         if(lakonl(1:1).eq.'F') cycle
         if(lakonl(1:7).eq.'DCOUP3D') cycle
!
         if(lakonl(7:8).ne.'LC') then
!
            imat=ielmat(1,i)
            amat=matname(imat)
            if(norien.gt.0) then
               iorien=ielorien(1,i)
            else
               iorien=0
            endif
!
            if(nelcon(1,imat).lt.0) then
               ihyper=1
            else
               ihyper=0
            endif
         else
!     
!     composite materials
!     
            mint2d=4
            nopes=8
!     determining the number of layers
!     
            nlayer=0
            do k=1,mi(3)
               if(ielmat(k,i).ne.0) then
                  nlayer=nlayer+1
               endif
            enddo
!     
!     determining the layer thickness and global thickness
!     at the shell integration points
!     
            iflag=1
            indexe=ipkon(i)
            do kk=1,mint2d
               xi=gauss3d2(1,kk)
               et=gauss3d2(2,kk)
               call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag)
               tlayer(kk)=0.d0
               do k=1,nlayer
                  thickness=0.d0
                  do j=1,nopes
                     thickness=thickness+thicke(k,indexe+j)*shp2(4,j)
                  enddo
                  tlayer(kk)=tlayer(kk)+thickness
                  xlayer(k,kk)=thickness
               enddo
            enddo
            iflag=3
!     
            ilayer=0
            do k=1,4
               dlayer(k)=0.d0
            enddo
!     
         endif
!     
         indexe=ipkon(i)
c     Bernhardi start
         if(lakonl(1:5).eq.'C3D8I') then
            nope=11
         elseif(lakonl(4:5).eq.'20') then
c     Bernhardi end
            nope=20
         elseif(lakonl(4:4).eq.'2') then
            nope=26
         elseif(lakonl(4:4).eq.'8') then
            nope=8
         elseif(lakonl(4:5).eq.'10') then
            nope=10
         elseif(lakonl(4:5).eq.'14') then
            nope=14
         elseif(lakonl(4:4).eq.'4') then
            nope=4
         elseif(lakonl(4:5).eq.'15') then
            nope=15
         elseif(lakonl(4:4).eq.'6') then
            nope=6
         elseif((lakonl(1:1).eq.'E').and.(lakonl(7:7).ne.'F')) then
!
!           spring elements, contact spring elements and
!           dashpot elements
!
            if(lakonl(7:7).eq.'C') then
!
!              contact spring elements
!
               if(mortar.eq.1) then
!
!                 face-to-face penalty
!
                  nope=kon(ipkon(i))
               elseif(mortar.eq.0) then
!
!                 node-to-face penalty
!
                  nope=ichar(lakonl(8:8))-47
                  konl(nope+1)=kon(indexe+nope+1)
               endif
            else
!
!              genuine spring elements and dashpot elements
!
               nope=ichar(lakonl(8:8))-47
            endif
         else
            cycle
         endif
!
         if(lakonl(4:5).eq.'8R') then
            mint3d=1
         elseif(lakonl(4:7).eq.'20RB') then
            if((lakonl(8:8).eq.'R').or.(lakonl(8:8).eq.'C')) then
               mint3d=50
            else
               call beamintscheme(lakonl,mint3d,ielprop(i),prop,
     &              null,xi,et,ze,weight)
            endif
         elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'26R').or.
     &          (lakonl(4:6).eq.'20R')) then
            if(lakonl(7:8).eq.'LC') then
               mint3d=8*nlayer
            else
               mint3d=8
            endif
         elseif(lakonl(4:4).eq.'2') then
            mint3d=27
         elseif((lakonl(4:5).eq.'10').or.(lakonl(4:5).eq.'14')) then
            mint3d=4
         elseif(lakonl(4:4).eq.'4') then
            mint3d=1
         elseif(lakonl(4:5).eq.'15') then
            mint3d=9
         elseif(lakonl(4:4).eq.'6') then
            mint3d=2
         elseif(lakonl(1:1).eq.'E') then
            mint3d=0
         endif

         do j=1,nope
            konl(j)=kon(indexe+j)
            do k=1,3
               xl(k,j)=co(k,konl(j))
               vl(k,j)=vold(k,konl(j))
               voldl(k,j)=vold(k,konl(j))
            enddo
         enddo
!
!     computation of the objective function and the derivate 
!
       if(desvar.gt.0) then 
          do ii=1,nope
             if(konl(ii).eq.nodedesi(desvar)) then
                xl(ndirdesi(desvar),ii)=xl(ndirdesi(desvar),ii)+distmin
                exit
             endif
          enddo
       endif
!
!        q contains the nodal forces per element; initialisation of q
!
         if((iperturb(1).ge.2).or.((iperturb(1).le.0).and.(iout.lt.1))) 
     &      then
            do m1=1,nope
               do m2=0,mi(2)
                  q(m2,m1)=fn(m2,konl(m1))
               enddo
            enddo
         endif
!
!        calculating the forces for the contact elements
!
         if(mint3d.eq.0) then
!
!           "normal" spring and dashpot elements
!
            kode=nelcon(1,imat)
            if(lakonl(7:7).eq.'A') then
               t0l=0.d0
               t1l=0.d0
               if(ithermal(1).eq.1) then
                  t0l=(t0(konl(1))+t0(konl(2)))/2.d0
                  t1l=(t1(konl(1))+t1(konl(2)))/2.d0
               elseif(ithermal(1).ge.2) then
                  t0l=(t0(konl(1))+t0(konl(2)))/2.d0
                  t1l=(vold(0,konl(1))+vold(0,konl(2)))/2.d0
               endif
            endif
!
!           spring elements (including contact springs)
!     
            if(lakonl(2:2).eq.'S') then
               if((lakonl(7:7).eq.'A').or.((mortar.eq.0).and.
     &          ((nmethod.ne.1).or.(iperturb(1).ge.2).or.(iout.ne.-1))))
     &               then
                  call springforc_n2f(xl,konl,vl,imat,elcon,nelcon,elas,
     &              fnl,ncmat_,ntmat_,nope,lakonl,t1l,kode,elconloc,
     &              plicon,nplicon,npmat_,ener(1,i),iener,
     &              stx(1,1,i),mi,springarea(1,konl(nope+1)),nmethod,
     &              ne0,nstate_,xstateini,xstate,reltime,
     &              ielas,ener(1,i+ne))
               elseif((mortar.eq.1).and.
     &           ((nmethod.ne.1).or.(iperturb(1).ge.2).or.(iout.ne.-1)))
     &               then
                  iloc=kon(indexe+nope+1)
                  jfaces=kon(indexe+nope+2)
                  igauss=kon(indexe+nope+1)
                  call springforc_f2f(xl,vl,imat,elcon,nelcon,elas,
     &              fnl,ncmat_,ntmat_,nope,lakonl,t1l,kode,elconloc,
     &              plicon,nplicon,npmat_,ener(1,i),iener,
     &              stx(1,1,i),mi,springarea(1,iloc),nmethod,
     &              ne0,nstate_,xstateini,xstate,reltime,
     &              ielas,iloc,jfaces,igauss,pslavsurf,pmastsurf,
     &              clearini,ener(1,i+ne))
               endif
!             
!              next lines are not executed in linstatic.c before the
!              setup of the stiffness matrix (i.e. nmethod=1 and
!              iperturb(1)<1 and iout=-1).
!
               if((lakonl(7:7).eq.'A').or.
     &           ((nmethod.ne.1).or.(iperturb(1).ge.2).or.(iout.ne.-1)))
     &              then
                  do j=1,nope
                     do k=1,3
                        fn(k,konl(j))=fn(k,konl(j))+fnl(k,j)
                     enddo
                  enddo
               endif
!
!              dashpot elements (including contact dashpots)
!
            elseif((nmethod.eq.4).or.(nmethod.eq.5).or.
     &             ((abs(nmethod).eq.1).and.(iperturb(1).ge.2))) then
            endif
         elseif(ikin.eq.1) then
            do j=1,nope
               do k=1,3
                  veoldl(k,j)=veold(k,konl(j))
               enddo
            enddo            
         endif
!
         do jj=1,mint3d
!
!        ---------------------------------------------------------------
!        Passing of stresses and strains from previous step 
!        for sensitivity analysis
!        ---------------------------------------------------------------
!
            do m1=1,6
               stiini(m1,jj,i)=sti(m1,jj,i)
            enddo
 
            if(lakonl(4:5).eq.'8R') then
               xi=gauss3d1(1,jj)
               et=gauss3d1(2,jj)
               ze=gauss3d1(3,jj)
               weight=weight3d1(jj)
            elseif(lakonl(4:7).eq.'20RB') then
               if((lakonl(8:8).eq.'R').or.(lakonl(8:8).eq.'C')) then
                  xi=gauss3d13(1,jj)
                  et=gauss3d13(2,jj)
                  ze=gauss3d13(3,jj)
                  weight=weight3d13(jj)
               else
                  call beamintscheme(lakonl,mint3d,ielprop(i),prop,
     &                 jj,xi,et,ze,weight)
               endif
            elseif((lakonl(4:4).eq.'8').or.
     &             (lakonl(4:6).eq.'20R').or.(lakonl(4:6).eq.'26R'))
     &        then
               if(lakonl(7:8).ne.'LC') then
                  xi=gauss3d2(1,jj)
                  et=gauss3d2(2,jj)
                  ze=gauss3d2(3,jj)
                  weight=weight3d2(jj)
               else
                  kl=mod(jj,8)
                  if(kl.eq.0) kl=8
!
                  xi=gauss3d2(1,kl)
                  et=gauss3d2(2,kl)
                  ze=gauss3d2(3,kl)
                  weight=weight3d2(kl)
!
                  ki=mod(jj,4)
                  if(ki.eq.0) ki=4
!
                  if(kl.eq.1) then
                     ilayer=ilayer+1
                     if(ilayer.gt.1) then
                        do k=1,4
                           dlayer(k)=dlayer(k)+xlayer(ilayer-1,k)
                        enddo
                     endif
                  endif
                  ze=2.d0*(dlayer(ki)+(ze+1.d0)/2.d0*xlayer(ilayer,ki))/
     &                 tlayer(ki)-1.d0
                  weight=weight*xlayer(ilayer,ki)/tlayer(ki)
!
!                 material and orientation
!
                  imat=ielmat(ilayer,i)
                  amat=matname(imat)
                  if(norien.gt.0) then
                     iorien=ielorien(ilayer,i)
                  else
                     iorien=0
                  endif
!     
                  if(nelcon(1,imat).lt.0) then
                     ihyper=1
                  else
                     ihyper=0
                  endif
               endif
            elseif(lakonl(4:4).eq.'2') then
               xi=gauss3d3(1,jj)
               et=gauss3d3(2,jj)
               ze=gauss3d3(3,jj)
               weight=weight3d3(jj)
            elseif((lakonl(4:5).eq.'10').or.(lakonl(4:5).eq.'14')) then
               xi=gauss3d5(1,jj)
               et=gauss3d5(2,jj)
               ze=gauss3d5(3,jj)
               weight=weight3d5(jj)
            elseif(lakonl(4:4).eq.'4') then
               xi=gauss3d4(1,jj)
               et=gauss3d4(2,jj)
               ze=gauss3d4(3,jj)
               weight=weight3d4(jj)
            elseif(lakonl(4:5).eq.'15') then
               xi=gauss3d8(1,jj)
               et=gauss3d8(2,jj)
               ze=gauss3d8(3,jj)
               weight=weight3d8(jj)
            elseif(lakonl(4:4).eq.'6') then
               xi=gauss3d7(1,jj)
               et=gauss3d7(2,jj)
               ze=gauss3d7(3,jj)
               weight=weight3d7(jj)
            endif
!
c     Bernhardi start
            if(lakonl(1:5).eq.'C3D8R') then
               call shape8hr(xl,xsj,shp,gs,a)
            elseif(lakonl(1:5).eq.'C3D8I') then
               call shape8hu(xi,et,ze,xl,xsj,shp,iflag)
            elseif(nope.eq.20) then
c     Bernhardi end
               if(lakonl(7:7).eq.'A') then
                  call shape20h_ax(xi,et,ze,xl,xsj,shp,iflag)
               elseif((lakonl(7:7).eq.'E').or.
     &                (lakonl(7:7).eq.'S')) then
                  call shape20h_pl(xi,et,ze,xl,xsj,shp,iflag)
               else
                  call shape20h(xi,et,ze,xl,xsj,shp,iflag)
               endif
            elseif(nope.eq.26) then
               call shape26h(xi,et,ze,xl,xsj,shp,iflag,konl)
            elseif(nope.eq.8) then
               call shape8h(xi,et,ze,xl,xsj,shp,iflag)
            elseif(nope.eq.10) then
               call shape10tet(xi,et,ze,xl,xsj,shp,iflag)
            elseif(nope.eq.14) then
               call shape14tet(xi,et,ze,xl,xsj,shp,iflag,konl)
            elseif(nope.eq.4) then
               call shape4tet(xi,et,ze,xl,xsj,shp,iflag)
            elseif(nope.eq.15) then
               call shape15w(xi,et,ze,xl,xsj,shp,iflag)
            else
               call shape6w(xi,et,ze,xl,xsj,shp,iflag)
            endif
!
!                 vkl(m2,m3) contains the derivative of the m2-
!                 component of the displacement with respect to
!                 direction m3
!
            do m2=1,3
               do m3=1,3
                  vkl(m2,m3)=0.d0
               enddo
            enddo
!
            do m1=1,nope
               do m2=1,3
                  do m3=1,3
                     vkl(m2,m3)=vkl(m2,m3)+shp(m3,m1)*vl(m2,m1)
                  enddo
c                  write(*,*) 'vnoeie',ivkl,konl(m1),(vkl(m2,k),k=1,3)
               enddo
            enddo
!
!           for frequency analysis or buckling with preload the
!           strains are calculated with respect to the deformed
!           configuration
!           for a linear iteration within a nonlinear increment:
!           the tangent matrix is calculated at strain at the end
!           of the previous increment
!
            if((iperturb(1).eq.1).or.(iperturb(1).eq.-1))then
               do m2=1,3
                  do m3=1,3
                     vokl(m2,m3)=0.d0
                  enddo
               enddo
!
               do m1=1,nope
                  do m2=1,3
                     do m3=1,3
                        vokl(m2,m3)=vokl(m2,m3)+
     &                       shp(m3,m1)*voldl(m2,m1)
                     enddo
                  enddo
               enddo
            endif
!
            kode=nelcon(1,imat)
!
!           calculating the strain
!
!           attention! exy,exz and eyz are engineering strains!
!
            exx=vkl(1,1)
            eyy=vkl(2,2)
            ezz=vkl(3,3)
            exy=vkl(1,2)+vkl(2,1)
            exz=vkl(1,3)+vkl(3,1)
            eyz=vkl(2,3)+vkl(3,2)
!
            if(iperturb(2).eq.1) then
!     
!                 Lagrangian strain
!     
               exx=exx+(vkl(1,1)**2+vkl(2,1)**2+vkl(3,1)**2)/2.d0
               eyy=eyy+(vkl(1,2)**2+vkl(2,2)**2+vkl(3,2)**2)/2.d0
               ezz=ezz+(vkl(1,3)**2+vkl(2,3)**2+vkl(3,3)**2)/2.d0
               exy=exy+vkl(1,1)*vkl(1,2)+vkl(2,1)*vkl(2,2)+
     &              vkl(3,1)*vkl(3,2)
               exz=exz+vkl(1,1)*vkl(1,3)+vkl(2,1)*vkl(2,3)+
     &              vkl(3,1)*vkl(3,3)
               eyz=eyz+vkl(1,2)*vkl(1,3)+vkl(2,2)*vkl(2,3)+
     &              vkl(3,2)*vkl(3,3)
!
!           for frequency analysis or buckling with preload the
!           strains are calculated with respect to the deformed
!           configuration
!
            elseif(iperturb(1).eq.1) then
               exx=exx+vokl(1,1)*vkl(1,1)+vokl(2,1)*vkl(2,1)+
     &              vokl(3,1)*vkl(3,1)
               eyy=eyy+vokl(1,2)*vkl(1,2)+vokl(2,2)*vkl(2,2)+
     &              vokl(3,2)*vkl(3,2)
               ezz=ezz+vokl(1,3)*vkl(1,3)+vokl(2,3)*vkl(2,3)+
     &              vokl(3,3)*vkl(3,3)
               exy=exy+vokl(1,1)*vkl(1,2)+vokl(1,2)*vkl(1,1)+
     &              vokl(2,1)*vkl(2,2)+vokl(2,2)*vkl(2,1)+
     &              vokl(3,1)*vkl(3,2)+vokl(3,2)*vkl(3,1)
               exz=exz+vokl(1,1)*vkl(1,3)+vokl(1,3)*vkl(1,1)+
     &              vokl(2,1)*vkl(2,3)+vokl(2,3)*vkl(2,1)+
     &              vokl(3,1)*vkl(3,3)+vokl(3,3)*vkl(3,1)
               eyz=eyz+vokl(1,2)*vkl(1,3)+vokl(1,3)*vkl(1,2)+
     &              vokl(2,2)*vkl(2,3)+vokl(2,3)*vkl(2,2)+
     &              vokl(3,2)*vkl(3,3)+vokl(3,3)*vkl(3,2)
            endif
!
c            if(iperturb(2).eq.1) then
c!     
c!                 Lagrangian strain
c!     
c               exx=exx+(vkl(1,1)**2+vkl(2,1)**2+vkl(3,1)**2)/2.d0
c               eyy=eyy+(vkl(1,2)**2+vkl(2,2)**2+vkl(3,2)**2)/2.d0
c               ezz=ezz+(vkl(1,3)**2+vkl(2,3)**2+vkl(3,3)**2)/2.d0
c               exy=exy+vkl(1,1)*vkl(1,2)+vkl(2,1)*vkl(2,2)+
c     &              vkl(3,1)*vkl(3,2)
c               exz=exz+vkl(1,1)*vkl(1,3)+vkl(2,1)*vkl(2,3)+
c     &              vkl(3,1)*vkl(3,3)
c               eyz=eyz+vkl(1,2)*vkl(1,3)+vkl(2,2)*vkl(2,3)+
c     &              vkl(3,2)*vkl(3,3)
c!     
c            endif
!
!              storing the local strains
!
            if(iperturb(1).ne.-1) then
               eloc(1)=exx
               eloc(2)=eyy
               eloc(3)=ezz
               eloc(4)=exy/2.d0
               eloc(5)=exz/2.d0
               eloc(6)=eyz/2.d0
            else
!
!              linear iteration within a nonlinear increment:
!
               eloc(1)=vokl(1,1)+
     &           (vokl(1,1)**2+vokl(2,1)**2+vokl(3,1)**2)/2.d0
               eloc(2)=vokl(2,2)+
     &           (vokl(1,2)**2+vokl(2,2)**2+vokl(3,2)**2)/2.d0
               eloc(3)=vokl(3,3)+
     &           (vokl(1,3)**2+vokl(2,3)**2+vokl(3,3)**2)/2.d0
               eloc(4)=(vokl(1,2)+vokl(2,1)+vokl(1,1)*vokl(1,2)+
     &              vokl(2,1)*vokl(2,2)+vokl(3,1)*vokl(3,2))/2.d0
               eloc(5)=(vokl(1,3)+vokl(3,1)+vokl(1,1)*vokl(1,3)+
     &              vokl(2,1)*vokl(2,3)+vokl(3,1)*vokl(3,3))/2.d0
               eloc(6)=(vokl(2,3)+vokl(3,2)+vokl(1,2)*vokl(1,3)+
     &              vokl(2,2)*vokl(2,3)+vokl(3,2)*vokl(3,3))/2.d0
            endif
!
!                 calculating the deformation gradient (needed to
!                 convert the element stiffness matrix from spatial
!                 coordinates to material coordinates
!                 deformation plasticity)
!
            if((kode.eq.-50).or.(kode.le.-100)) then
!
!                    calculating the deformation gradient
!
c     Bernhardi start
               xkl(1,1)=vkl(1,1)+1.0d0
               xkl(2,2)=vkl(2,2)+1.0d0
               xkl(3,3)=vkl(3,3)+1.0d0
c     Bernhardi end
               xkl(1,2)=vkl(1,2)
               xkl(1,3)=vkl(1,3)
               xkl(2,3)=vkl(2,3)
               xkl(2,1)=vkl(2,1)
               xkl(3,1)=vkl(3,1)
               xkl(3,2)=vkl(3,2)
!
!                    calculating the Jacobian
!
               vj=xkl(1,1)*(xkl(2,2)*xkl(3,3)-xkl(2,3)*xkl(3,2))
     &              -xkl(1,2)*(xkl(2,1)*xkl(3,3)-xkl(2,3)*xkl(3,1))
     &              +xkl(1,3)*(xkl(2,1)*xkl(3,2)-xkl(2,2)*xkl(3,1))
!
!              inversion of the deformation gradient (only for
!              deformation plasticity)
!
               if(kode.eq.-50) then
!
                  ckl(1,1)=(xkl(2,2)*xkl(3,3)-xkl(2,3)*xkl(3,2))/vj
                  ckl(2,2)=(xkl(1,1)*xkl(3,3)-xkl(1,3)*xkl(3,1))/vj
                  ckl(3,3)=(xkl(1,1)*xkl(2,2)-xkl(1,2)*xkl(2,1))/vj
                  ckl(1,2)=(xkl(1,3)*xkl(3,2)-xkl(1,2)*xkl(3,3))/vj
                  ckl(1,3)=(xkl(1,2)*xkl(2,3)-xkl(2,2)*xkl(1,3))/vj
                  ckl(2,3)=(xkl(2,1)*xkl(1,3)-xkl(1,1)*xkl(2,3))/vj
                  ckl(2,1)=(xkl(3,1)*xkl(2,3)-xkl(2,1)*xkl(3,3))/vj
                  ckl(3,1)=(xkl(2,1)*xkl(3,2)-xkl(2,2)*xkl(3,1))/vj
                  ckl(3,2)=(xkl(3,1)*xkl(1,2)-xkl(1,1)*xkl(3,2))/vj
!
!                 converting the Lagrangian strain into Eulerian
!                 strain (only for deformation plasticity)
!
                  cauchy=0
                  call str2mat(eloc,ckl,vj,cauchy)
               endif
!                        
            endif
!
!                 calculating fields for incremental plasticity
!
            if(kode.le.-100) then
!
!              calculating the deformation gradient at the
!              start of the increment
!
!              calculating the displacement gradient at the
!              start of the increment
!
               do m2=1,3
                  do m3=1,3
                     vikl(m2,m3)=0.d0
                  enddo
               enddo
!
               do m1=1,nope
                  do m2=1,3
                     do m3=1,3
                        vikl(m2,m3)=vikl(m2,m3)
     &                       +shp(m3,m1)*vini(m2,konl(m1))
                     enddo
                  enddo
               enddo
!
!              calculating the deformation gradient of the old
!              fields
!
               xikl(1,1)=vikl(1,1)+1
               xikl(2,2)=vikl(2,2)+1.
               xikl(3,3)=vikl(3,3)+1.
               xikl(1,2)=vikl(1,2)
               xikl(1,3)=vikl(1,3)
               xikl(2,3)=vikl(2,3)
               xikl(2,1)=vikl(2,1)
               xikl(3,1)=vikl(3,1)
               xikl(3,2)=vikl(3,2)
!
!              calculating the Jacobian
!
               vij=xikl(1,1)*(xikl(2,2)*xikl(3,3)
     &              -xikl(2,3)*xikl(3,2))
     &              -xikl(1,2)*(xikl(2,1)*xikl(3,3)
     &              -xikl(2,3)*xikl(3,1))
     &              +xikl(1,3)*(xikl(2,1)*xikl(3,2)
     &              -xikl(2,2)*xikl(3,1))
!
!              stresses at the start of the increment
!
               do m1=1,6
                  stre(m1)=stiini(m1,jj,i)
               enddo
!
            endif
!
!                 prestress values
!
            if(iprestr.eq.1) then
               do kk=1,6
                  beta(kk)=-prestr(kk,jj,i)
               enddo
            else
               do kk=1,6
                  beta(kk)=0.d0
               enddo
            endif
!
            if(ithermal(1).ge.1) then
!
!              calculating the temperature difference in
!              the integration point
!
               t0l=0.d0
               t1l=0.d0
               if(ithermal(1).eq.1) then
                  if((lakonl(4:5).eq.'8 ').or.
     &               (lakonl(4:5).eq.'8I')) then
                     do i1=1,8
                        t0l=t0l+t0(konl(i1))/8.d0
                        t1l=t1l+t1(konl(i1))/8.d0
                     enddo
                  elseif((lakonl(4:6).eq.'20 ').or.
     &                   (lakonl(4:6).eq.'26 ')) then
                     nopered=20
                     call lintemp(t0,t1,konl,nopered,jj,t0l,t1l)
                  else
                     do i1=1,nope
                        t0l=t0l+shp(4,i1)*t0(konl(i1))
                        t1l=t1l+shp(4,i1)*t1(konl(i1))
                     enddo
                  endif
               elseif(ithermal(1).ge.2) then
                  if((lakonl(4:5).eq.'8 ').or.
     &               (lakonl(4:5).eq.'8I')) then
                     do i1=1,8
                        t0l=t0l+t0(konl(i1))/8.d0
                        t1l=t1l+vold(0,konl(i1))/8.d0
                     enddo
                  elseif((lakonl(4:6).eq.'20 ').or.
     &                   (lakonl(4:6).eq.'26 ')) then
                     nopered=20
                     call lintemp_th(t0,vold,konl,nopered,jj,t0l,t1l,mi)
                  else
                     do i1=1,nope
                        t0l=t0l+shp(4,i1)*t0(konl(i1))
                        t1l=t1l+shp(4,i1)*vold(0,konl(i1))
                     enddo
                  endif
               endif
               tt=t1l-t0l
            endif
!
!                 calculating the coordinates of the integration point
!                 for material orientation purposes (for cylindrical
!                 coordinate systems)
!
            if((iorien.gt.0).or.(kode.le.-100)) then
               do j=1,3
                  pgauss(j)=0.d0
                  do i1=1,nope
                     pgauss(j)=pgauss(j)+shp(4,i1)*co(j,konl(i1))
                  enddo
               enddo
            endif
!
!                 material data; for linear elastic materials
!                 this includes the calculation of the stiffness
!                 matrix
!
            istiff=0
!     
            call materialdata_me(elcon,nelcon,rhcon,nrhcon,alcon,
     &           nalcon,imat,amat,iorien,pgauss,orab,ntmat_,
     &           elas,rho,i,ithermal,alzero,mattyp,t0l,t1l,ihyper,
     &           istiff,elconloc,eth,kode,plicon,nplicon,
     &           plkcon,nplkcon,npmat_,plconloc,mi(1),dtime,i,jj,
     &           xstiff,ncmat_)
!
!           determining the mechanical strain
!
            if(ithermal(1).ne.0) then
               do m1=1,6
                  emec(m1)=eloc(m1)-eth(m1)
c                 emec0(m1)=emeini(m1,jj,i)
               enddo
            else
               do m1=1,6
                  emec(m1)=eloc(m1)
c                 emec0(m1)=emeini(m1,jj,i)
               enddo
            endif
            if(kode.le.-100) then
               do m1=1,6
                  emec0(m1)=emeini(m1,jj,i)
               enddo
            endif
!
!           subtracting the plastic initial strains
!
            if(iprestr.eq.2) then
               do m1=1,6
                  emec(m1)=emec(m1)-prestr(m1,jj,i)
               enddo
            endif
!
!           calculating the local stiffness and stress
!
            call mechmodel(elconloc,elas,emec,kode,emec0,ithermal,
     &           icmd,beta,stre,xkl,ckl,vj,xikl,vij,
     &           plconloc,xstate,xstateini,ielas,
     &           amat,t1l,dtime,time,ttime,i,jj,nstate_,mi(1),
     &           iorien,pgauss,orab,eloc,mattyp,qa(3),istep,iinc,
     &           ipkon,nmethod,iperturb)
!
            if(((nmethod.ne.4).or.(iperturb(1).ne.0)).and.
     &         (nmethod.ne.5)) then
               do m1=1,21
                  xstiff(m1,jj,i)=elas(m1)
               enddo
            endif
!
            if(iperturb(1).eq.-1) then
!
!                    if the forced displacements were changed at
!                    the start of a nonlinear step, the nodal
!                    forces due do this displacements are 
!                    calculated in a purely linear way, and
!                    the first iteration is purely linear in order
!                    to allow the displacements to redistribute
!                    in a quasi-static way (only applies to
!                    quasi-static analyses (*STATIC))
!
               eloc(1)=exx-vokl(1,1)
               eloc(2)=eyy-vokl(2,2)
               eloc(3)=ezz-vokl(3,3)
               eloc(4)=exy-(vokl(1,2)+vokl(2,1))
               eloc(5)=exz-(vokl(1,3)+vokl(3,1))
               eloc(6)=eyz-(vokl(2,3)+vokl(3,2))
!
               if(mattyp.eq.1) then
                  e=elas(1)
                  un=elas(2)
                  um=e/(1.d0+un)
                  al=un*um/(1.d0-2.d0*un)
                  um=um/2.d0
                  am1=al*(eloc(1)+eloc(2)+eloc(3))
                  stre(1)=am1+2.d0*um*eloc(1)
                  stre(2)=am1+2.d0*um*eloc(2)
                  stre(3)=am1+2.d0*um*eloc(3)
                  stre(4)=um*eloc(4)
                  stre(5)=um*eloc(5)
                  stre(6)=um*eloc(6)
               elseif(mattyp.eq.2) then
                  stre(1)=eloc(1)*elas(1)+eloc(2)*elas(2)
     &                 +eloc(3)*elas(4)
                  stre(2)=eloc(1)*elas(2)+eloc(2)*elas(3)
     &                 +eloc(3)*elas(5)
                  stre(3)=eloc(1)*elas(4)+eloc(2)*elas(5)
     &                 +eloc(3)*elas(6)
                  stre(4)=eloc(4)*elas(7)
                  stre(5)=eloc(5)*elas(8)
                  stre(6)=eloc(6)*elas(9)
               elseif(mattyp.eq.3) then
                  stre(1)=eloc(1)*elas(1)+eloc(2)*elas(2)+
     &                 eloc(3)*elas(4)+eloc(4)*elas(7)+
     &                 eloc(5)*elas(11)+eloc(6)*elas(16)
                  stre(2)=eloc(1)*elas(2)+eloc(2)*elas(3)+
     &                 eloc(3)*elas(5)+eloc(4)*elas(8)+
     &                 eloc(5)*elas(12)+eloc(6)*elas(17)
                  stre(3)=eloc(1)*elas(4)+eloc(2)*elas(5)+
     &                 eloc(3)*elas(6)+eloc(4)*elas(9)+
     &                 eloc(5)*elas(13)+eloc(6)*elas(18)
                  stre(4)=eloc(1)*elas(7)+eloc(2)*elas(8)+
     &                 eloc(3)*elas(9)+eloc(4)*elas(10)+
     &                 eloc(5)*elas(14)+eloc(6)*elas(19)
                  stre(5)=eloc(1)*elas(11)+eloc(2)*elas(12)+
     &                 eloc(3)*elas(13)+eloc(4)*elas(14)+
     &                 eloc(5)*elas(15)+eloc(6)*elas(20)
                  stre(6)=eloc(1)*elas(16)+eloc(2)*elas(17)+
     &                 eloc(3)*elas(18)+eloc(4)*elas(19)+
     &                 eloc(5)*elas(20)+eloc(6)*elas(21)
               endif
            endif
! 
!           calculation of the shape energy for the objective 
!           function: SHAPE ENERGY
!
            ener(jj,i)=enerini(jj,i)+
     &            ((eloc(1)-eth(1)-emeini(1,jj,i))*
     &        (stre(1)+stiini(1,jj,i))+
     &        (eloc(2)-eth(2)-emeini(2,jj,i))*
     &        (stre(2)+stiini(2,jj,i))+
     &        (eloc(3)-eth(3)-emeini(3,jj,i))*
     &        (stre(3)+stiini(3,jj,i)))/2.d0+
     &         (eloc(4)-eth(4)-emeini(4,jj,i))*(stre(4)+stiini(4,jj,i))+
     &         (eloc(5)-eth(5)-emeini(5,jj,i))*(stre(5)+stiini(5,jj,i))+
     &         (eloc(6)-eth(6)-emeini(6,jj,i))*(stre(6)+stiini(6,jj,i))

            enerint=enerint+weight*xsj*ener(jj,i)
!
!        end of loop over all integration points
         enddo
!
!        q contains the contributions to the nodal force in the nodes
!        belonging to the element at stake from other elements (elements
!        already treated). These contributions have to be
!        subtracted to get the contributions attributable to the element
!        at stake only
!
         if(calcul_qa.eq.1) then
            do m1=1,nope
               do m2=1,3
                  qa(1)=qa(1)+dabs(fn(m2,konl(m1))-q(m2,m1))
               enddo
            enddo
            nal=nal+3*nope
         endif
! 
!     summation of the objective function over all elements
!
      dg=dg+enerint
! 
!     end of loop over all elements in thread
!
      enddo
!
!        for the nominal model ( without variation of the nodes) stop after 
!        sum over all elements once
!
         if(desvar.eq.0) then
            g0(numobject)=dg
         else
!     
!        finite difference step 
!
            dgdx(desvar,numobject)=(dg-g0(numobject))/distmin         
         endif
!
!     end of loop over all designvariables
!
      enddo
!     
      return
      end
